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PACS. 83.60.Rs - Shear rate-dependent structure (shear thinning and shear thickening). 

PACS. 83.60.Wc - Flow instabilities. 

PACS. 82.40.Bj - Oscillations, chaos, and bifurcations. 



Abstract. - We study the strain response to steady imposed stress in a spatially homogeneous, 
scalar model for shear thickening, in which the local rate of yielding P(Z) of mesoscopic 'elastic 
elements' is not monotonic in the local strain I. Despite this, the macroscopic, steady-state flow 
curve (stress vs. strain rate) is monotonic. However, for a broad class of P(/), the response to 
steady stress is not in fact steady flow, but spontaneous oscillation. We discuss this finding in 
relation to other theoretical and experimental flow instabilities. Within the parameter ranges 
we studied, the model does not exhibit rheo-chaos. 



The flow behaviour of shear-thickening materials such as dense colloidal suspensions can be 
complex For example, imposition of a steady mean strain rate can lead to large, possibly 

chaotic, variations in the mean stress pj] . The same occurs in some types of shear-thickening 
micellar surfactant solutions, where true temporal chaos seems now to be established || (and 
also in shear thinning systems; see 0). Other unexpected behaviour, such as a bifurcation to 
an oscillatory state, has also been seen in shear-thickening 'onion' phases of surfactant ||. It is 
not yet known to what extent such unsteady flow is generic in shear-thickening systems; in this 
letter we attempt to shed some light on the issue by studying a much-simplified, generic model. 
In this model we find, for a wide range of parameters, spontaneous rheological oscillation of 
the strain rate at fixed stress. Rheo-chaos is, however, not found for the parameters studied 
so far. 

A feature that distinguishes the rheological instabilities encountered in shear-thickening 
from those arising in Newtonian fluids is that the nonlinearity is not inertial (not from the 
advective term of the Navier Stokes equation): the Reynolds number is essentially zero [0]. 
Instead it arises from anharmonic elastic responses at large deformations, complicated and 
perhaps strongly enhanced by the presence of memory effects. Flow instabilities leading 
to chaos have been studied recently by Grosso et al. in a model for suspended rodlike 
particles. As that work shows, and our work confirms, temporal instabilities can arise even in 
a model where macroscopic spatial inhomogeneity is disallowed altogether. This is a strong 
demarcation from the familiar shear-banding instabilities that arise whenever the steady-state 
flow curve is nonmonotonic (the flow curve is the function 17(7), where a is shear stress and 7 
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the rate of shear strain) j^, ^) . Shear banding is well documented in shear thinning materials 
(particularly micelles || ) but also possible in shear thickening ones , albeit with a different 
spatial organization of the bands of coexisting material. However it is not the subject of this 
Letter - we are interested in temporal inhomogeneity. 

Our simplified model of a shear thickening fluid is defined as follows. We imagine an 
ensemble of mesoscopic elastic elements, each having a local strain variable I. We consider 
only simple shear strains, and neglect normal stresses, so that the only nontrivial stress is 
the corresponding shear stress a. (These assumptions reduce the problem to a scalar one.) 
Let us define P(l,t) to be the probability density function of elements which have a local 
strain I at time t. We assume P(l,t) evolves in time according to two distinct mechanisms: 
homogeneous shearing at a rate 7, and the yielding of elements at a rate T(l) per unit time. 
In calculating the stress, the elements are supposed to behave elastically between yield events, 
so that the local stress is just kl, where k > is an elastic constant. The global stress a 
is simply the arithmetic mean of the local stresses, or a = (kl) = k(l). (Here the angled 
brackets represent an instantaneous average over P(l, t).) We ignore completely the fact that, 
in practice, the strain rate 7 can vary locally in space in response to each mesoscopic element 
having a different local stress. 

In the present work, T(l) is the same for all elements, which simplifies the analysis signif- 
icantly. A more comprehensive model, in which elements each have their own 'yield strain' 
E is studied in a longer paper fll|] , where further details of some of our calculations can also 
be found. Both models are closely related to the (shear-thinning) SGR model of Sollich et 
al. H^l, but we will choose a very different form of T(l) from that of SGR. 

The master equation for P(l,t) is then 

DP dP 

^+i-Qj- = -r(i)p + u,(t)6(i) (1) 

The second term on the left-hand side represents the increase in local strains / according 
to the spatially uniform strain rate 7. The two terms on the right hand side describe the 
birth and death of elements with strain I, respectively, where 5(1) reflects the assumption that 
newly-yielded elements are unstrained. The total rate of yielding u>(t) is defined by 

/oo 
dlT(l)P(l,t) = (T(l)) (2) 
-00 

The key remaining ingredient of our model is the hopping rate T(l). We parameterise this 
as follows: 

= T exp [- (E - kl 2 /2) /x(l)) (3) 

This is a pseudo-activated form, in which each element is subject to an effective temperature 
x, and attempts (at a rate To) to hop over a barrier of height E ~ ^kl 2 . The latter expression 
allows for the lowering of the elastic yield barrier by the imposition of strain; on its own this 
will always lead to shear thinning. 

However, the novel feature in the current model is to allow x(l) to be a decreasing function 
of /. This reflects the intuition that, even if a small local stress or strain always promotes 
yield, a large enough one may jam an element against its neighbours causing the jump rate to 
fall. Other than in steady state (see below) we have been unable to find an analytic solution 
to (|l|) for any non-trivial x(l). Instead it has been numerically integrated. With the above as 
motivation, we chose for numerical purposes relatively simple (piecewise constant) decreasing 
functions x(l). The precise choice of x(l) — in particular, whether it is smooth or not — does 
not seem to qualitatively influence the results. 
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Results. - The steady state solution Poo (I) of (|l]) is found as (setting k = T = 1 from 
now on, for convenience) 

PooQ) = UooT 1 exp [-7 _1 /(0] (4) 

where f(l) = J^dlT(l'), and Woo, the asymptotic jump rate, is fixed by normalization of Poo- 
It is straightforward to show from (Q) that in the limit of slow flows 7 — > + , the steady 
state stress response Coo = (I) 00 ( m an obvious notation) is always linear: (Too ~ ^e E / xty0 \ 
Hence there is no yield stress for any choice of x(0). This contrasts with a model having an 
exponential distribution of barrier heights E for different elements, which does show onset 
of a yield stress, connected with the presence of a glass transition, as x is reduced |Q. For 
monodisperse E, as here, there is no such transition. 

We now show that the steady state flow curve, for any choice of function x(l), has a 
monotonically increasing (7(7). First, by differentiating (l}oo and using (|j), we find 

da 

OO OO Coo 1 ,,„\ /_n 

-^r- = ^ ; I" To ('//oo ( 5 ) 

97 Woo 07 7 7 
Similarly, the normalisation integral for Poo can be differentiated with respect to 7 to give 

1 dujoc 1 1 /A 

w-:^7" = ^-^ (/> - (6) 

Combining these two expressions (with <7oo = (Ooo) produces 

= ('/>ao-(0oo(/)oc (7) 

which is positive since /(Z) is monotone increasing. Thus there are no regions on the flow curve 
with negative slope, and for no choice of x(l) does one expect any shear-banding instability 
to arise. However, the argument does not prove that a steady state is ever actually achieved. 
We have found that it is reached under conditions of imposed strain j(t), but not under an 
imposed stress a (for appropiate x(l)), for which the model exhibits temporal oscillations 
instead. 

This numerical finding (detailed below) is also supported by a linear stability analysis for 
perturbations about the steady state flow, described fully in |ll]] . Such an analysis shows that 
the transient behaviour close to steady state is never a real exponential decay but always has 
an oscillatory component. When the decay rate of these oscillations changes sign, the flow 
becomes unstable to permanent oscillation. 

The oscillatory regime. - The time evolution of the system at fixed stress was studied 
numerically by methods described in JTlJ] . (The problem is not standard and significant 
amounts of computer time were required to achieve reliable results.) A range of different 
functions x(l) and of E values were tried. Varying a over a wide range of values in each case 
reveals that steady flow is always reached for sufficiently small and sufficiently large imposed 
stresses; but for intermediate a, an oscillatory regime is observed whenever x(l) is strongly 
enough decreasing. It is impractical to explore the entire function space for x(l) so at present 
the precise criteria for the presence of oscillations are not known. 

Within the oscillatory regime, we have so far found only single-period oscillations with no 
sign of period doubling or other more complex nonsteady behaviour that might point to the 
onset of rheo-chaos in some parts of the parameter space. Some examples of the oscillatory 
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behaviour are given in Fig. |T](a). To obtain these results we took x = 1 for I < 0.4 and x = 0.4 
for I > 0.4; in addition we set E — 5. The strain j(t) is shown for three different imposed a. 
On close inspection the oscillations, after a brief transient, appear perfectly periodic. This 
applies not only in the stress but in the underlying distribution P(l,t), which determines the 
future evolution. Hence there seems little doubt that these oscillations will persist indefinitely. 

From Fig. |l](a) it is remarkable that the mean strain rate 7, defined as 7 averaged over 
a single period of oscillation, is clearly a decreasing function of a, in complete contrast to 
the monotonic, steady state, flow curve. The time averaged values 7 as read off from the 
simulations are plotted against a in Fig. 0(b), overlayed with the steady state flow curve. 
Within the oscillatory regime, the 7 line deviates from the flow curve to the extent that it 
loses monotonicity over a large interval of a. If the oscillations were not detected but only the 
average behaviour measured, this would resemble the flow curve of a discontinuously shear- 
thickening system close to a jamming transition, although in that case shear banding could 
be expected as well. (Such a flow curve is indeed found in models where x depends not on 
the local stress I but only on its global mean, a; see flT|| .) 

Upper and lower transition points between oscillatory and steady flow can be seen in the 
figure. At each transition point a = a c , the period of oscillation remains finite whilst the 
amplitude vanishes smoothly according to \a — cr c \ a with a > 0. As we discuss in our longer 
work JTl| , technical difficulties have prevented us from reliably fixing the asymptotic values 
of a. In particular we have been unable to rule out a = 0.5 for both transitions, as expected 
for a Hopf bifurcation [^|,|l3| (although since the control parameter is an integral constraint, 
it is not clear to us if it should be a Hopf bifurcation). 

Though always perfectly periodic, the waveform of the oscillations varies with a. Close to 
either transition, the oscillations are near-sinusoidal, as demonstrated in Fig. |^(a). Further 
into the oscillatory regime, 7 can no longer be decomposed into a single harmonic, but instead 
approaches a waveform in which most of the variation in 7 is compressed into a small fraction 
of the total period of oscillation. An example is given in Fig. ||(b) . This behaviour at first 
appears to be similar to the 'stick-slip' motion observed in systems such as earthquakes ]l5j, 



ultra-thin liquid films |l6|] and granular media 17, 18[. However, the underlying physics in 
our model seems to be somewhat different to these examples. Indeed, stick-slip is usually 
viewed as a surface phenomenon, whereas the model studied in this paper has no surface and 
describes only bulk effects. Irrespective of the waveform, the product of 7 and the period of 
oscillation T is approximately constant, I* — 7T « 2.3 for this example. This will be explained 
below in terms of the evolution of P(l, t). 

Mechanism of oscillation. Snapshots of P{l,t) during a single period of oscillation 
are given in Fig. [|. The mechanism behind the oscillations can be qualitatively described 
in terms of two coexisting populations of elements — a 'hot' population of unstrained or 
slightly strained elements with I ks 0, and hence a high effective temperature x(l); and a 'cold' 
population of elements with I 3> and therefore a low x(l). Starting at a time when j(t) 
is near to is minimum value, corresponding to the first snapshot in the figure, the stress- 
weighted yield rate {IT (I)) for both populations of elements is low. This is because the hot 
elements have small I, whereas the cold elements have a low yield rate T{1). Thus the strain 
rate 7 required to maintain the imposed constant mean stress a, which is = {IT (I)) (as 
seen by multiplying (Q) by Idl and integrating), is also low. 

Although 7(f) is small in this state, it is nonetheless non-zero, and therefore the strains I 
of the cold elements increase according to I = 7. This decreases their effective energy barrier 
E — 1 2 /2, increasing their hop rate and therefore the rate of stress lost due to to cold elements 
yielding. Thus j(t) = {IP {I)) will also increase, which accelerates the rate at which the cold 
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elements yield, and so on. This description is that of a positive feedback loop that causes 7 
to increase at ever faster rates until all of the cold elements have been depleted. The second 
snapshot in Fig. || show the state of the system shortly before this has happened. 

While j(t) is high, the hot elements are becoming rapidly strained and have little time 
to yield before crossing over into the cold region with large I. Thus P(l,t) is flat for small /. 
As 7(£) again decreases, this flat part of the distribution will start to decay as the elements 
within it yield. However, the yield rate depends on x(l), and since this changes to a lower 
value at I = 0.4, P(l,t) will decay more rapidly for I smaller than 0.4 than for I greater than 
0.4. Thus a dip will occur around the point I w 0.4, which can clearly be seen in Fig. ||. This 
dip becomes more pronounced with time until the system can again be viewed as coexisting 
'hot' and 'cold' populations, when the cycle begins again. 

A possible interpretation of the parameter /*, defined previously, is that it is the amount 
by which the system needs to be strained until the positive feedback loop just described starts 
to dominate the system behaviour, causing it to 'reset' to the start of its cycle. If this is 
correct, then /* should correspond to the point at which highly-strained elements have the 
same yield rate as unstrained elements, i.e. T(0) « T(7*). Rough calculations based on this 
assumption give I* « y/2E[l — x(oo)/a;(0)], which for this example predicts /* w 2.4, in fair 
agreement with the observed value. 

Discussion. - Within a simple model of shear thickening, we have found regimes where 
steady flow at constant imposed stress is unstable and temporal oscillation occurs instead. 
Known instances of such rheological instability have often been explained in terms of the 
spatial coexistence of subpopulations, or phases. For instance, the temporal oscillations in 
viscosity observed in wormlike micelles under an imposed stress was attributed to a slowly 
fluctuating interface between a fluid phase and shear-induced structures pj . For surfactant 
solutions in the lyotropic lamellar phase, it was attributed to coexisting ordered and disordered 
phases [pi. By contrast, the temporal oscillations observed in our model arise even though they 
are by assumption spatially homogeneous, and the flow curves are everywhere monotonic. It 
may be that, in practice, such rheo-oscillations will always be coupled to spatial heterogeneity 
of some sort. Experiments are not advanced enough to indicate whether this is the case. 

One distinguishing feature of our equations is the role of memory effects, rather than 
inertia, in allowing oscillations. In our model, the memory resides in P(l,t) which can have 
different shapes for the same applied stress (I). There may thus be a link to work on 'delay 
differential equations' in which a first order differential equation for a state variable contains a 
feedback term that depends on the same state variable at an earlier time [jl9| . Such equations 
are capable of showing a range of instabilities, including chaos, despite the absence of inertia. 
The nature of this link remains a topic of current research. Its exploration might help answer 
the following question: what features, if any, could be added to our simple rheological model of 
shear thickening to allow rheo-chaos to arise without coupling to spatial degrees of freedom |?J . 
This could in turn help answer another important question: what is the dimension of the 
strange attractors that arise in rheo-chaos ||? In conclusion we suggest that the study of 
temporal (or spatiotemporal) pattern formation in non-Newtonian flows (at effectively zero 
Reynolds number) will form a major topic of experimental and theoretical research in the 
coming years. 
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Fig. 1 - (a) The strain y(t) under an imposed stress a = 0.1 (solid line), 0.13 (dashed line) and 0.2 
(dot-dashed line) for x{V) — 1 for I < 0.4 and x — 0.4 for I > 0.4. The system was initially unstrained 
and E = 5. (b) The mean strain rate j against the imposed stress a for the same system. The 
solid and open circles represent steady and oscillatory solutions, respectively, as observed from the 
simulations. The sizes of the circles are larger than the error bars. For comparison, the theoretical 
flow curve for the steady state solution M) is plotted as a solid line. 




Fig. 2 - Examples of the oscillatory strain response to a constant stress for the same x(l) as in Fig. |l| 
The imposed stresses are (a) a = 0.075, which is just above the threshold value for steady flow, and 
(b) a — 0.4, which is well into the oscillatory regime. Note that the axes in (b) are semilogarithmic. 
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Fig. 3 - Snapshots of P(l, i) in a monodisperse system with E\ — 5, and x(l) — 1 for I < 0.4, 0.4 
otherwise. The imposed stress was a = 0.2. The times of each plot are t = 1.1 x 10 5 (solid line), 
t — 1.25 x 10 5 (dashed line) and t — 1.3 x 10 5 (dot-dashed line), corresponding to before, during and 
just after the peak in 7, respectively. The vertical dotted line is where x(l) changes value. Note that 
this corresponds to the lower line in Fig. hi. 



